Task7_Redkokosh
Redkokosh Kirill
12/7/2021
library(factoextra)
library(cluster)diabetes <- read.csv("winequality-red.csv", header = TRUE, as.is = FALSE, sep = ';')
diabetes <- diabetes[, -1]
head(diabetes)## class glucose insulin sspg
## 1 Normal 80 356 124
## 2 Normal 97 289 117
## 3 Normal 105 319 143
## 4 Normal 90 356 199
## 5 Normal 90 323 240
## 6 Normal 86 381 157Первичный анализ данных
Описательная статистика
summary(diabetes)## class glucose insulin sspg
## Chemical:36 Min. : 70 Min. : 45.0 Min. : 10.0
## Normal :76 1st Qu.: 90 1st Qu.: 352.0 1st Qu.:118.0
## Overt :33 Median : 97 Median : 403.0 Median :156.0
## Mean :122 Mean : 540.8 Mean :186.1
## 3rd Qu.:112 3rd Qu.: 558.0 3rd Qu.:221.0
## Max. :353 Max. :1568.0 Max. :748.0Среднее и медиана отличается у всех признаков, ожидаем несимметричные распределения.
Виды признаков
Построим matrix plot для того, чтобы увидеть особенности в наших данных.
library(lattice)
library(ggplot2)
library('GGally')
ggpairs(diabetes, title="correlogram", columns=c(2:4), upper = list(continuous = "points"), diag = list(continuous = "barDiag"))
Прологарифмируем признаки.
diabetes_l <- transform(diabetes, glucose=log(glucose), insulin=log(insulin))
names(diabetes_l)[names(diabetes_l) == 'glucose'] <- 'log_glucose'
names(diabetes_l)[names(diabetes_l) == 'insulin'] <- 'log_insulin'
ggpairs(diabetes_l, title="correlogram", columns=c(2:4), upper = list(continuous = "points"), diag = list(continuous = "barDiag"))
Удалим единичные outliers.
А также добавим раскраску по признаку “class”.
diabetes_lo <- diabetes_l
diabetes_lo[rownames(diabetes_lo)[diabetes_lo$log_insulin < 5],] <- NA
diabetes_lo <- na.omit(diabetes_lo)
ggpairs(diabetes_lo, title="correlogram", columns=c(2:4), upper = list(continuous = "points"), diag = list(continuous = "barDiag"), mapping=ggplot2::aes(colour = class))
При дальнейшей кластеризации классифицирующий признак– “class” рассматриваться не будет.
library("FactoMineR")
dflo <- diabetes_lo[, -1]
res.pca <- PCA(dflo, scale.unit = TRUE, ncp = 2, graph= FALSE)
get_eigenvalue(res.pca)## eigenvalue variance.percent cumulative.variance.percent
## Dim.1 2.10928604 70.309535 70.30953
## Dim.2 0.83660334 27.886778 98.19631
## Dim.3 0.05411062 1.803687 100.00000fviz_pca_biplot(res.pca, habillage=diabetes_lo$class)
Первые 2 ГК описывают 98% данных.
Рассмотрим 3 метода для определения числа кластеров:
- «Силуэт» вычисляется на основе среднего внутрикластерного расстояния и среднего расстояния до ближайшего кластера по каждому образцу.
Для точки \(i\) из кластера \(C_{I}\):
\(a(i)={\frac{1}{|C_{I}|-1}}\sum _{j\in C_{I},i\neq j}d(i,j)\) – среднее расстояние между i-той точкой и всеми остальными точками кластера.
\(b(i)=\min _{J\neq I}{\frac {1}{|C_{J}|}}\sum _{j\in C_{J}}d(i,j)\) – наименьшее среднее расстояние между i-той точкой и всеми точками любого другого кластера.
Теперь мы определяем силуэт (значение) одной точки данных:
\(s(i) = \begin{cases} 1-a(i)/b(i), & \mbox{if } a(i) < b(i) \\ 0, & \mbox{if } a(i) = b(i) \\ b(i)/a(i)-1, & \mbox{if } a(i) > b(i) \\ \end{cases}\)
Среднее s(i) по всем точками кластера - это мера того, насколько плотно сгруппированы все точки кластера. Таким образом, среднее s(i) по всем данным всего набора данных является мерой того, насколько правильно были кластеризованы данные.
- Метод локтя основывается на сумме квадратов расстояний между точками и центроидами.
Рассматривается характер изменения общего внутригруппового разброса с увеличением числа групп k. Объединив все n наблюдений в одну группу, мы имеем наибольшую внутрикластерную дисперсию, которая будет снижаться до 0 при \(k \rightarrow n\). На каком-то этапе снижение этой дисперсии замедляется - на графике это происходит в точке, называемой “локтем”.
- Статистика разрыва сравнивает общую дисперсию внутри кластера для различных значений k с их ожидаемыми значениями для распределения без кластеризации.
Пусть \(E^∗_n{log(W^∗_k)}\) обозначает оценку средней дисперсии \(W^∗_k\), полученной бутстреп-методом, когда k кластеров образованы случайными наборами объектов из исходной выборки размером n.
Тогда статистика:
\(Gap_n(k)=E^∗_n{log(W^∗_k)}−log(W_k)\)
определяет отклонение наблюдаемой дисперсии \(W_k\) от ее ожидаемой величины при справедливости нулевой гипотезы о том, что исходные данные образуют только один кластер.
При сравнительном анализе последовательности значений \(Gap_n(k),k=2,…,K_{max}\) наибольшее значение статистики соответствует наиболее полезной группировке, дисперсия которой максимально меньше внутригрупповой дисперсии кластеров, собранных из случайных объектов исходной выборки.
fviz_nbclust(dflo, kmeans, method = "silhouette", k.max = 10) #"silhouette" (for average silhouette width), "wss" (for total within sum of square) and "gap_stat" (for gap statistics).
fviz_nbclust(dflo, kmeans, method = "wss", k.max = 10)
fviz_nbclust(dflo, kmeans, method = "gap_stat", k.max = 6)
Будем рассматривать деление на 2 и 3 класстерa.
set.seed(10)
km2 <- kmeans(dflo, centers = 2, nstart = 100)
km2$size## [1] 117 27km2$centers#/sapply(dflo, sd)## log_glucose log_insulin sspg
## 1 4.74357 6.185327 139.9145
## 2 4.60564 6.145251 378.7037km3 <- kmeans(dflo, centers = 3, nstart = 100)
km3$size## [1] 52 9 83km3$centers#/sapply(dflo, sd)## log_glucose log_insulin sspg
## 1 4.562684 6.006067 239.4615
## 2 4.640043 6.213764 534.6667
## 3 4.823253 6.281513 112.4217Кластеры получаются неравного размера (как и при рассмотрении классифицирующего признака, но там другое соотвношение индивидов в классах).
Деление на кластеры в плоскости первых 2 ГК и на ggpairs для 2 и 3 кластеров при использовании метода k-means:
fviz_cluster(km2, data = dflo)
fviz_cluster(km3, data = dflo)
ggpairs(dflo, title="correlogram", columns=c(1:3), upper = list(continuous = "points"), diag = list(continuous = "barDiag"), mapping=ggplot2::aes(colour = factor(km2$cluster)))
ggpairs(dflo, title="correlogram", columns=c(1:3), upper = list(continuous = "points"), diag = list(continuous = "barDiag"), mapping=ggplot2::aes(colour = factor(km3$cluster)))
Кластеризируем данные с меньшим числом признаков, для этого рассмотрим первые 2 ГК (на одну меньше, чем количество кластеров).
km3_pca <- kmeans(res.pca$ind$coord, centers = 3, nstart = 100)
km3_pca$size## [1] 18 26 100km3_pca$centers#/sapply(res.pca$ind$coord, sd)## Dim.1 Dim.2
## 1 -0.8814062 1.8586132
## 2 2.8245783 0.1323704
## 3 -0.5757372 -0.3689667fviz_cluster(km3_pca, data = res.pca$ind$coord)
#ggpairs(data.frame(res.pca$ind$coord), title="correlogram", columns=c(1:2), upper = list(continuous = "points"), diag = list(continuous = "barDiag"), mapping=ggplot2::aes(colour = factor(km3_pca$cluster)))Деление также неравномерное.
Перейдем к следующему методу кластеризации– методу иерархической кластеризации с разными правилами объединения кластеров.
library(dplyr)
d <- dist(scale(dflo), method = "euclidean")
h_average <- hclust(d, method = "average")
plot(h_average, cex = 0.7, hang = -1) #агломеративная кластеризация, групповое среднее расстояние
cut_average <- cutree(h_average, k = 3)
df_cl <- mutate(dflo, cluster = cut_average)
ggpairs(df_cl, title="correlogram", columns=c(1:3), upper = list(continuous = "points"), diag = list(continuous = "barDiag"), mapping=ggplot2::aes(colour = factor(cluster)))
h_single <- hclust(d, method = "single")
plot(h_single, cex = 0.7, hang = -1) #агломеративная кластеризация, расстояние ближнего соседа (цепочки)
cut_single <- cutree(h_single, k = 3)
df_cl <- mutate(dflo, cluster = cut_single)
ggpairs(df_cl, title="correlogram", columns=c(1:3), upper = list(continuous = "points"), diag = list(continuous = "barDiag"), mapping=ggplot2::aes(colour = factor(cluster)))
Тут все совсем плохо ):
res.hc <- hclust(d, method = "complete") #агломеративная кластеризация, расстояние дальнего соседа (шарики)
plot(res.hc, cex = 0.7, hang = -1)
rect.hclust(res.hc, k = 3, border = 2:5)
cut_complete <- cutree(res.hc, k = 3)
df_cl <- mutate(dflo, cluster = cut_complete)
ggpairs(df_cl, title="correlogram", columns=c(1:3), upper = list(continuous = "points"), diag = list(continuous = "barDiag"), mapping=ggplot2::aes(colour = factor(cluster)))
Дивизионная кластеризация diana:
На шаге 0 все объекты объединены в один кластер.
На каждом шаге кластер делится, пока на шаге n - 1 все объекты данных не будут отделены друг от друга.
На каждом шаге делится кластер, назовем его R на два кластера A и B. Первоначально A равно R и B пуст. На первом этапе мы должны переместить один объект из A в B. Для каждого объекта i из A мы вычисляем среднее несходство ко всем другим объектам A:
\(d(i,A({i}))=\frac{1}{|A|-1}\sum_{i \neq j}d(i,j)\)
Объект i’, для которого уравнение выше достигает своего максимального значения, будет перемещено, поэтому положим
\(A_{new}=A_{old} \setminus \{i'\}\)
\(B_{new}=B_{old} \bigcup\{i'\}\)
На следующих этапах мы ищем другие точки для перехода из A в B. Пока A все еще содержит более одного объекта, мы вычисляем
\(d(i,A({i}))-d(i,B)=\frac{1}{|A|-1}\sum_{j \in A, i \neq j}d(i,j) - \frac{1}{|B|}\sum_{h \in B}d(i,h)\)
для каждого объекта i из A, и мы рассматриваем объект i’‘, который максимизирует эту величину. Когда максимальное значение вышеприведенного уравнения строго положительное, мы перемещаем i’’ от A до B, а затем просматриваем \(A_{new}\). С другой стороны, когда максимальное значение разности отрицательное или 0, мы останавливаем процесс, и разделение R на A и B завершено.
На каждом шаге делящего алгоритма мы также должны решить, какой кластер нужно разбить. Для этого вычислим диаметр
\(diam(Q)=max_{j \in Q, h \in Q} d(j,h)\)
для каждого кластера Q, доступного после предыдущего шага, и выбираем кластер, с наибольшим диаметром.
divisive.clust <- diana(as.matrix(d), diss = TRUE, keep.diss = TRUE) #дивизионная кластеризация
plot(divisive.clust, main = "Divisive", which.plots = 2, hang = -1)
cut_div <- cutree(divisive.clust, k = 3)
df_cl <- mutate(dflo, cluster = cut_div)
ggpairs(df_cl, title="correlogram", columns=c(1:3), upper = list(continuous = "points"), diag = list(continuous = "barDiag"), mapping=ggplot2::aes(colour = factor(cluster)))
Многие методы показывают похожие результаты. Если сравнивать с исходной классификацией, то методы кластеризации разбивают данные на другие группы.
divisive.clust <- diana(as.matrix(dist(scale(res.pca$ind$coord))), diss = TRUE, keep.diss = TRUE)
plot(divisive.clust, main = "Divisive", which.plots = 2, hang = -1)
cut_div <- cutree(divisive.clust, k = 3)
df_cl <- mutate(dflo, cluster = cut_div)
ggpairs(df_cl, title="correlogram", columns=c(1:3), upper = list(continuous = "points"), diag = list(continuous = "barDiag"), mapping=ggplot2::aes(colour = factor(cluster)))
Не слишком хорошо сработал метод понижения числа признаков.
library(mclust)
mc <- mclust::Mclust(dflo)
summary(mc) #первый символ относится к объему, второй - к форме, а третий - к ориентации. (E-- равный; V-- переменный; I-- расположение относительно осей координат)## ----------------------------------------------------
## Gaussian finite mixture model fitted by EM algorithm
## ----------------------------------------------------
##
## Mclust VVV (ellipsoidal, varying volume, shape, and orientation) model with 3
## components:
##
## log-likelihood n df BIC ICL
## -701.7836 144 29 -1547.692 -1559.786
##
## Clustering table:
## 1 2 3
## 33 83 28head(mc$z) #каждому наблюдению назначается кластер с максимальной оцененной вероятностью.## [,1] [,2] [,3]
## 1 4.352532e-03 0.9956475 1.496784e-14
## 2 1.344389e-07 0.9999999 8.739538e-39
## 3 8.770285e-07 0.9999991 4.940825e-32
## 4 2.443657e-03 0.9975563 2.713471e-15
## 5 5.391036e-04 0.9994609 2.013616e-19
## 6 9.513377e-03 0.9904866 7.771781e-12par(mfrow = c(1, 2))
plot(mc, "classification")
plot(mc, "uncertainty")
fviz_cluster(mc, data = dflo)
Результаты похожи на исходную классификацию.
Максимальный BIC (среди всех моделей) у модели VVV (но у неё же и максимальное число параметров).
Попробуем выбрать другую модель:
plot(mc, "BIC") #Зависимость критерия BIC от числа кластеров для различных вариантов параметризации ковариационной матрицы
apply(mc$BIC, 2, which.max)## EII VII EEI VEI EVI VVI EEE VEE EVE VVE EEV VEV EVV VVV
## 8 9 6 5 4 4 7 6 4 4 3 3 3 3Кроме первых двух моделей BIC у остальных моделей близок, возьмем модель, кластеризующую данные на 3 кластера, но имеющюю меньше параметров.
mclustModelNames('EVV')## $model
## [1] "EVV"
##
## $type
## [1] "ellipsoidal, equal volume"mc1 <- mclust::Mclust(dflo, modelNames = 'EVV')
summary(mc1)## ----------------------------------------------------
## Gaussian finite mixture model fitted by EM algorithm
## ----------------------------------------------------
##
## Mclust EVV (ellipsoidal, equal volume) model with 3 components:
##
## log-likelihood n df BIC ICL
## -707.3937 144 27 -1548.972 -1560.778
##
## Clustering table:
## 1 2 3
## 25 88 31head(mc1$z) ## [,1] [,2] [,3]
## 1 1.830318e-04 0.9998170 1.069044e-15
## 2 9.980166e-16 1.0000000 1.527052e-41
## 3 9.212325e-14 1.0000000 1.761841e-34
## 4 6.984441e-05 0.9999302 2.128655e-16
## 5 2.713809e-07 0.9999997 7.528674e-21
## 6 1.465955e-03 0.9985340 8.978124e-13par(mfrow = c(1, 2))
plot(mc1, "classification")
plot(mc1, "uncertainty")
fviz_cluster(mc1, data = dflo)
Результаты немного отличаются, но все еще похожи на исходную классификацию.